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Abstract. The localized deposition of the energy of a laser pulse, as it ablates a 
solid target, introduces high thermal pressure gradients in the plasma. The thermal 
expansion of this laser-heated plasma into the ambient medium (ionized residual gas) 
triggers the formation of non-linear structures in the collision-less plasma. Here 
an electron-proton plasma is modelled with a particle-in-cell (PIC) simulation to 
reproduce aspects of this plasma expansion. A jump is introduced in the thermal 
pressure of the plasma, across which the otherwise spatially uniform temperature and 
density change by the factor 100. The electrons from the hot plasma expand into the 
cool one and the charge imbalance drags a beam of cool electrons into the hot plasma. 
This double layer reduces the electron temperature gradient. The presence of the low- 
pressure plasma modifies the proton dynamics compared to the plasma expansion into 
a vacuum. The jump in the thermal pressure develops into a primary shock. The fast 
protons, which move from the hot into the cold plasma in form of a beam, give rise 
to the formation of phase space holes in the electron and proton distributions. The 
proton phase space holes develop into a secondary shock that thermalizes the beam. 
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1. Introduction 

The impact of a laser pulse on a solid target results in the evaporation of the target 
material. The heated plasma expands under its own thermal pressure and shocks as 
well as other nonlinear plasma structures form. Generating collision-less plasma shocks 
in a laboratory experiment permits us to study their detailed dynamics in a controlled 
manner. A better understanding of such shocks is not only relevant for the laser- 
plasma experiment as such and for inertial confinement fusion experiments. It can also 
provide further insight into the dynamics of solar system shocks and the nonrelativistic 
astrophysical shocks, like the supernova remnant shocks [H El El HI [5] . 

An obstacle to an in-depth investigation of the laser-generated shocks has been 
so far, that the frequently used optical probing techniques could not resolve the shock 
structure at the required spatio-temporal resolution. The now available proton imaging 
technique [61 [7] helps us overcoming this limitation. This method can provide accurate 
spatial electric field profiles at a high time resolution, as long as no strong magnetic 
fields are present. The nonrelativistic flow speed of the laser-generated shock, e.g. that 
in Ref. [8j, implies that no strong self- induced magnetic fields due to the filamentation 
instability or the mixed mode instability [9], \TU\ occur at the shock front. 

The availability of electric field data at a high resolution serves as a motivation 
to perform related numerical simulations and to compare their results with the 
experimental ones. The experimental observations from Ref. [8] , which are most relevant 
for the simulation study we perform here, can be summarized as follows. The ablation 
of a solid target consisting of aluminium or tungsten by a laser pulse with a duration of 
~ 470 ps and an intensity of 10 15 W/cm 2 results in a plasma with a density w 10 18 cm -3 
and with an electron temperature of a few keV. This plasma expands into an ambient 
plasma with the density < 10 15 cm~ 3 . The ambient plasma has been produced mainly 
by a photo-ionization of the residual gas. The dominant components of the residual 
gas, which consists of diluted air, are oxygen and nitrogen. Electrostatic structures, 
which move through the ionized residual g observed. Their propagation speeds 

suggest that one is an electrostatic shock [11] with a thickness of a few electron Debye 
lengths, which is expanding approximately with the ion acoustic velocity 2 — 4 x 10 5 m/s. 
Ion-acoustic solitons are trailing the shock. Another structure moves at twice the shock 
speed, which is probably related to a shock-reflected ion beam. The electron-electron, 
electron-ion and ion-ion mean free paths for the residual gas have been determined for 
this particular experiment. They are of the order of cm and much larger than the shock 
width of a few tens of /im. The shock and the electrostatic structures are collision-less. 

The experiment can measure the electric fields, the propagation speed of the electric 
field structures and it can estimate the electron temperature and density. The bulk 
parameters of the ions, such as their temperature, mean speed and ionization state, 
are currently inaccessible, as well as detailed information about the spatial distribution 
of the plasma. We can set up a plasma simulation with the experimentally known 
parameters, and we can introduce an idealized model for the unknown initial conditions. 
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The detailed information about the state of the plasma, which is provided by Vlasov 
simulations [12J or by particle-in-cell (PIC) simulations [T3|,[TI], can then provide further 
insight into the expansion of this plasma. 

Here we investigate a mechanism that could result in the shock observed in Ref. [8]. 
We model with PIC simulations the interplay of two plasmas with a large difference in 
the thermal pressure, which are initially spatially separated. We aim at determining the 
spatio-temporal scale, over which a shock forms under this initial assumption, and we 
want to reveal the structures that develop in the wake of the shock. The temperature and 
density of the hot laser-ablated plasma both exceed initially that of the cold ambient 
plasma by two orders of magnitude. The density ratio is less than that between the 
expanding and the ambient plasma in Ref. [Sj. However, the density will not change 
in form of a single jump in the experiment and realistic density changes will probably 
be less or equal to the one we employ. Selecting the same jump in the density and 
temperature is computationally efficient, because both plasmas have the same Debye 
length that determines the grid cell size and the allowed time step. The ion temperature 
in the experiment is likely to be less than that of the electrons. The electron distribution 
can also not be approximated by two separate spatially uniform and thermal electron 
clouds, because the plasma generation is not fast compared to the electron diffusion. 
We show, however, that the shock forms long after the electrons have diffused in the 
simulation box and reached almost the same temperature everywhere. 

A change in the thermal pressure by a factor 10 4 should imply a plasma expansion 
that is similar to that into a vacuum. This process has received attention in the context 
of auroral, astrophysical and laser-generated plasmas and it has been investigated 
analytically within the framework of fluid models [T5J, [16] or Vlasov models [171 f!8] . 
It has been modelled numerically using a cold ion fluid and Boltzmann-distributed 
electrons [19] and with kinetic Vlasov and PIC simulations [201 l2l]. The plasma 
expansion of hot electrons and cool ions into a tenuous medium has also been examined 
with PIC simulations, such as the pioneering study in Ref. [22], which reported the 
formation of a double layer [231 EH [25] that cannot form if the plasma expands into a 
vacuum. Our simulation examines also the dynamics of protons as a first step towards a 
simulation of the mix of oxygen and nitrogen ions that constitute the residual gas in the 
physical experiment. Notable differences between the expansion of the hot and dense 
plasma into the ambient plasma and the expansion into a vacuum are observed. 

The structure of this paper is the following. We describe the PIC method in Section 
2 and we give the initial conditions and the simulation parameters. Section 3 models 
the initial phase of the plasma expansion at a high phase space resolution, revealing 
details of the electron expansion and of the quasi-equilibrium, which is established for 
the electrons. A double layer develops at the thermal pressure jump, which drags the 
electrons from the tenuous plasma into the hot plasma in form of a cool beam. The 
electrons from the hot plasma leak into the cool plasma, which reduces the temperature 
difference between both plasmas. Section 4 examines the proton dynamics. The ambient 
plasma modifies the proton expansion. The thermal pressure jump evolves into a shock, 
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which moves approximately with the proton thermal speed of the hot plasma. If the 
plasma expands into a vacuum, then a plasma density change can only be accomplished 
by ion beams [21J, while the plasma is here compressed by the shock. The fastest 
protons in our simulation form a beam that outruns the shock. It interacts with the 
protons of the ambient medium to form phase space holes in the electron and proton 
distributions. The proton phase space holes develop into a secondary shock ahead of 
the primary one. This process may result in secondary shocks in experiments, similar 
to the radiation-driven ones [26]. The results are summarized in Section 5. 

2. The PIC simulation method and the initial conditions 

A PIC code approximates a plasma by an ensemble of computational particles (CPs), 
each of which is representing a phase space volume element. Each CP follows a phase 
space trajectory that is determined through the Lorentz force equation by the electric 
field E(x, t) and the magnetic field B(x, t). Both fields are evolved self-consistently 
in time using the Maxwell's equations and the macroscopic current J(x, t), which is 
the sum over the microcurrents of all CPs. The standard PIC method considers only 
collective interactions between particles, although some collisional effects are introduced 
through the interaction of CPs with the field fluctuations [27] . 

Collision operators have been prescribed for PIC simulations [28l [29] . The 
structures in the addressed experiment form and evolve in a plasma, in which collisional 
effects are not strong and such operators are thus not introduced here. We may 
illustrate this with the help of the electron collision rate v e ~ 2.9 x 10~ 6 n e In A T~ 3 / 2 s _1 
and the ion collision rate v { rj 4.8 x lO -8 ^ 4 // -1 / 2 n. In AT^~ 3 ^ 2 s _1 [30J for a spatially 
uniform plasma with the number density n e = rti = 10 15 (cm" 3 ) and the temperature 
T e = Ti = 10 3 (eV). We take a Coulomb logarithm In A = 10 and we consider oxygen 
with = 16. Both collision rates are comparable, if the mean ion charge Z « 4. We 
assume v e « z/j. The electron plasma frequency u p ~ 10 12 s _1 gives the low relative 
collision frequency v e jio v ~ 10~ 6 . The plasma flow in the experiment and other aspects, 
which are not taken into account by this simplistic estimate, alter this collision frequency. 
The mean-free path has been estimated to be of the order of a cm [8] and the ion beam 
with the speed 4 x 10 5 m/s crosses this distance during the time u p t ~ 25000. This 
presumably forms the upper time limit, for which we can neglect collisions. 

The presence of particles with keV energies and the preferential expansion direction 
of the plasma in the experiment imply, that mult i- dimensional PIC simulations should 
be electromagnetic in order to resolve the potentially important magnetic Weibel 
instabilities, which are driven by thermal anisotropies [3T] . Such instabilities can grow 
in the absence of relativistic beams of charged particles, but they are typically weaker 
than the beam-driven ones [32J . Here we restrict our simulation to one spatial dimension 
x (ID) and we set B(x,t = 0). The plasma expands along x and all particle beams 
will have velocity vectors aligned with x. The magnetic beam-driven instabilities have 
wavevectors that are oriented obliquely or perpendicular to the beam velocity vector 
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and they are not resolved by a ID simulation. The wavevectors, which are destabilized 
by the Weibel instability, can be aligned with the simulation direction, but only if the 
plasma is cooler along x than orthogonally to it. Such a thermal anisotropy can probably 
not form. Our electromagnetic simulation confirms that no magnetic instability grows. 
The ratio of the magnetic to the total energy remains at noise levels below 10~ 4 . 

A ID PIC simulation should provide a reasonable approximation to those sections 
of the expanding plasma front observed in Ref. [8], which are planar over a sufficiently 
wide spatial interval orthogonal to the expansion direction. We set the length of the 
ID simulation box to L. The plasma 1 is consisting of electrons (species 1) and protons 
(species 2), each with the density rih and the temperature Th — 1 keV, and it fills up 
the half-space —L/2 < x < 0. A number density = 10 15 cm~ 3 should be appropriate 
with regard to the experiment. The half-space < x < L/2 is occupied by the plasma 2, 
which is composed of electrons (species 3) and protons (species 4) with the temperature 
T c = 10 eV and the density n c = n^/100. All plasma species have initially a Maxwellian 
velocity distribution, which is at rest in the simulation frame. 

The ablated target material drives the plasma expansion but its ions are probably 
not involved in the evolution of the shock and of the other plasma structures. These 
structures are observed already 100-200 ps after the laser impact at a distance of about 1 
mm from the target. Aluminium ions, which are with a mass tra the lightest constituents 
of the target material, would have the thermal speed (T/mA) 1 ^ 2 ~ 10 5 m/s for T = 1 
keV. Hundred times this speed or a temperature of 10 MeV would be necessary for them 
to propagate 1 mm in 0.1 ns. We thus assume here that the shock and the other plasma 
structures involve only the ions of the residual gas, which is air at a low pressure. If 
we assume that these ions have a high ionization state and comparable charge-to-mass 
ratios, then the protons may provide a reasonable approximation to their dynamics. 

The equations solved by the PIC code are normalized with the number density 
nh, the plasma frequency Qi = (n/ l e 2 /m e eo) 1 ^ 2 and the Debye length Ad = fti/Oi of 
species 1, which equals that of the other species. The thermal speeds of the respective 
species are v t j = (Tj/mj) 1 ^ 2 , where j is the species index. We express the charge q k 
and mass m k of the k th CP in units of the elementary charge e and electron mass m e . 
Quantities in physical units have the subscript p and we substitute E p = Qiv t im e 'E/e, 
B p = Oim e B/e, J p = ev a n h 3, p p = en h p, x p = X D x, t p = t/£li, v p = vv a and 
p p = m e m k v t ip. The ID PIC code solves with v n = v n /c the equations 

V x B = v 2 n (d t ~E + J) , V x E = -d t B, V • E = p, V ■ B = 0, (1) 
dtPk = qk (E[x fe ] + v fc x B[x k )) , d t x k = v kjX . (2) 

The Lorentz force is solved for each CP with index k, position x k and velocity v^. It is 
necessary to interpolate the electromagnetic fields from the grid to the particle position 
to update p k and the microcurrents of each CP have to be interpolated to the grid 
to update the electromagnetic fields. Interpolation schemes are detailed in Ref. [13J. 
Our code is based on the virtual particle electromagnetic particle-mesh method [14] and 
it uses the lowest possible interpolation order possible with this scheme. Our code is 
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parallelized through the distribution of the CPs over all processors. 

Simulation 1 (Section 3) resolves the box length Lg = 3350 by Ng = 5 x 10 3 grid 
cells of size A x g = 0.67Ad. The dense species 1 and 2 are each resolved by 8 x 10 4 
CPs per cell and the tenuous species 3 and 4 by 800 CPs per cell, respectively. The 
simulation is evolved in time for the duration tg = 800, subdivided into 45000 time steps 
Atg. Simulation 2 in Section 4 resolves the box length = 10 Lg by Nl = 2.5 x 10 4 grid 
cells of size A x l = 1.34Ad. This grid cell size is sufficiently small to avoid a significant 
numerical self-heating [33J of the plasma during the simulation time. The total energy 
in the simulation is preserved to within » 10~ 5 . The species 1 and 2 are approximated 
by 6400 CPs per cell each and the species 3,4 by 64 CPs per cell, respectively. The 
system is evolved during ti = 25500 with 6.4 x 10 5 time steps. 

We use periodic boundary conditions for the particles and the fields in all directions. 
Ideally, no particles or waves should traverse the full box length during the simulation 
duration. The group velocity for the electrostatic waves and the propagation speed of 
the electrons are both comparable to Vt\. We obtain vatg/Lg « 0.24 for simulation 1 
and Vf\ti,/Li, ~ 0.76 for simulation 2. Both simulations ran on 16 CPUs on an AMD 
Opteron cluster (2.2 GHz). Simulation 1/2 ran for 100/800 hours. 

3. Simulation 1: Initial development 

Our initial conditions involve a jump in the bulk plasma properties at x ~ 0. Some 
electrons of the plasma 1 will expand into the half-space x > occupied by the plasma 2. 
The slow protons can not keep up with the electrons and the resulting charge imbalance 
gives rise to an electrostatic field E x . This E x confines the electrons of plasma 1 and 
it accelerates the electrons from the plasma 2 into the half-space x < 0. The electrons 
of the plasma 1 and 2 with x < are separated along the velocity direction by the 
electrostatic potential and form a double layer. 

Figure [1] examines the E x and its potential. The amplitude of E x peaks initially 
at x ~ and it accelerates the electrons into the negative x-direction. The position of 
the maximum of E x moves to larger x with increasing times and the peak amplitude 
decreases. The spatial profile of E x is smooth, which contrasts the one that drives 
the plasma expansion into a vacuum that has a cusp [21J. The potential difference 
~ 5 kV between plasma 1 and 2 remains unchanged. The spatial interval, in which 
the amplitude of E x is well above noise levels, is bounded. An interesting property of 
the double layer can thus be inferred according to [25]. Its electrostatic field can only 
redistribute the momentum between the four plasma species, but it can not provide a 
net flow momentum. This is true if the double layer is one-dimensional and electrostatic. 
The decrease of the peak electric field in Fig. []](b) resembles that in Fig. 3 in Ref. [19|. 
The decreasing electric force, in turn, implies that the ion acceleration in the Fig. 4 of 
Ref. [19] decreases as the time progresses, which should hold for our simulation too. 

The plasma phase space distribution at t = 60 is investigated in Fig. [2j A tenuous 
hot beam of electrons is diffusing from the plasma 1 into the half-space x > 0, while the 
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Figure 1. The electric field: (a) shows E x at the times t = 60, 120, 180, 240 and 300. 
The maximum amplitude decreases with time as (b) is showing and the location of the 
electric field maximum moves towards positive x (c). The potential in kV obtained 
from the E x distributions from (a) is displayed in (d). The potential jump remains 
unchanged, but the gradient is eroded with the time. 




(b) Position (d) Position 



Figure 2. (Colour online) The plasma distribution at t = 60: (a) shows the electron 
phase space distribution. Most electrons from the dense plasma remain confined to 
x < 0, but some diffuse into the tenuous plasma, (b) shows the proton phase space 
distribution. Some protons with v x > are accelerated in < x < 5. The protons with 
x, v x < stream freely to lower values of x. (c) The electron phase space distribution 
reveals a double layer, (a-c) show the 10-logarithmic number of CPs. (d) shows the 
number of CPs per cell of the electrons (dashed curve) and the protons (solid curve). 
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mean speed of the electrons of the plasma 2 becomes negative. The electrons of plasma 
1 and 2 with x < are separated by a velocity gap of ~ v t i/10. The protons that were 
close to the initial boundary x = at t = have propagated until t = 60 for a distance, 
which is proportional to their speed. A sheared velocity distribution can thus be seen 
in Fig. [2(b). The fastest protons of the plasma 1 with x > have also been accelerated 
by the E x by about v t2 /2, reaching now a peak speed ~ 4t> t2 . The fastest protons are 
found to the right of the maximum of E x at x ~ 2 at t = 60 in FigQJa). A similar 
acceleration is observed for the protons of plasma 2 in < x < 5. The densities of the 
electrons and protons disagree in the interval —5 < x < 5 and the net charge results 
in the electrostatic field E x > 0. Both curves in Fig. [2(d) intersect at x ~ 2, which 
coincides with the position in Fig. Ufa), where the E x has its maximum at t — 60. 

The density of the cold protons in Ref. [21] is practically discontinuous at the front 
of the expanding plasma, while it changes smoothly in our simulation. This is a result 
of our high proton temperature, which causes the thermal diffusion of the protons. The 
contour lines of the electron phase space density are curved at x ~ 0. Most electrons of 
plasma 1 that move to increasing values of x are reflected by the electrostatic potential 
at x ~ 0. These density contour lines resemble those of the distribution of electrons 
that expand into a vacuum at an early time in Ref. [21], which are all reflected by the 
potential at the plasma front. Here the inflow of electrons from plasma 2 into plasma 
1 allows some of the electrons of plasma 1 to overcome the potential. The electrons 
provide all energy for the proton expansion in Ref. [21] and their distribution develops 
a flat top. Here the proton thermal energy is the main driver and consequently the 
electron velocity distribution shows no clear deviation from a Maxwellian at any time. 

Figure [3] shows the plasma phase space distributions at the times t = 120 and 
t = 180. The plasma distributions are qualitatively similar to that in Fig. [2j Electrons 
diffuse out from plasma 1 into plasma 2, forming a hot beam, while the electrons of 
plasma 2 are dragged into the half-space x < in form of a cold beam. The confined 
electrons of plasma 1 expand to increasing x at a speed, which is determined mainly 
by the protons. The proton distribution shows an increasing velocity shear, but the 
apparent phase space boundary between the protons of plasma 1 and 2 is still intersecting 
v x = at x = 0. The front of the protons of plasma 1 at t = 120 and t = 180 is close 
to the position of the maximum of E x in Fig. [T](a) at x ~ 5 for t = 120 and x ~ 10 for 
t = 180. The protons at the front of plasma 1 and the protons of plasma 2 in the same 
interval are accelerated by the E x > and reach the peak speed ~ 5vt2- 

The electrons of plasma 1 in Fig. 0] at t — 300 have expanded into the half-space 
x > for several hundred Debye lengths. The electrons from the plasma 2, which 
have been dragged towards x < 0, interact with the electrons of plasma 1 through a 
two-stream instability. A chain of large electron phase space holes has developed for 
—400 < x < —300, which thermalize the beam distribution. No two-stream instability 
is yet observed in the interval x > 0, even though a beam distribution is present, for 
example, at x ~ 250. The change of the mean speed of the electron beam leaked from 
plasma 1 for x > inhibits the resonance that gives rise to the two-stream instability. 
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(b) Position (d) Position 

Figure 3. (Colour online) The 10-logarithmic phase space densities in units of CPs: 
The electron distribution in (a) and the proton distribution in (b) are sampled at 
t = 120, while (c) and (d) show them at t — 180. The protons in the interval x, v x < 
convect almost freely away from x — 0. The protons of the dense plasma in x,v x > 
accelerate. Electrons diffuse from the plasma 1 into the plasma 2 and form a hot beam, 
while electrons from the plasma 2 enter the plasma 1 in form of a cold beam. 
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(b) Position " < d >V v ti 

Figure 4. (Colour online) The 10-logarithmic number of CPs representing the 
electrons (a) and the protons (b) at the time t = 300. The electrons of plasma 1 
have spread out to x w 700. The protons of plasma 1 with x > are accelerated to 
about 5vt2- The electron density in units of CPs for x < is displayed in (c). Electron 
phase space holes are present for x < —300. The electron distribution integrated over 
250 < x < 260 is shown in (d). 



The mean speed of the electrons of plasma 2 does not vanish any more and it varies 
along x > to provide the return current that cancels that of the electrons of plasma 
1. The E x has noticably accelerated the protons in the interval 10 < x < 30, which still 
show the sheared distribution in the interval —25 < x < 25. 

The evolution of the plasma is animated in the movies 1 (electrons) and 2 (protons). 
The axis labels v e h = v t \ and v p h = v t 2- The colour scale denotes the 10-logarithmic 
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Figure 5. (Colour online) The 10-logarithmic phase space distributions, normalized 
to their respective peak values: (a) shows the proton distribution in simulation 1 and 
(b) that in simulation 2. The electron distributions in the simulations 1 and 2 are 
displayed in (c) and (d), respectively. 



number of CPs. The movie 1 reveals that a thin band of electrons parallel to v x 
propagates away instantly from the plasma 1 and towards x > 0. These electrons 
leave the plasma 1, before the E x has grown. The electrons diffusing into x > at 
later times, when the E x has developed, form a tenuous beam with a broad velocity 
spread. The electrons of plasma 1 can overcome the double layer potential of ~ 5 kV 
if their speed is v > 3v t i prior to the encounter of its electrostatic field. The movie 
1 furthermore illustrates the growth of the two-stream instability between the electron 
beam originating from the plasma 2 and the confined electrons of plasma 1 in x < 
and its saturation through the formation of electron phase space holes. The movie 2 
demonstrates, how the velocity shear of the protons develops and how the fastest protons 
of plasma 1 in x > are accelerated by E x . Neither the Fig. H] nor the movie 2 reveal 
the formation of a shocked proton distribution prior to the time tg- 

We expand the simulation box and we reduce the statistical representation of the 
plasma. Ideally, the plasma evolution should be unchanged. Figure [5] compares the 
plasma data provided by simulation 1 (box length L s ) and by simulation 2 (L L = 10L S ) 
at the time t,s, when we stop simulation 1. The proton distributions in both simulations 
are practically identical and we notice only one quantitative difference. The sheared 
proton distribution of plasma 1 extends to x ~ —60 and v x ~ — 3v t 2 in simulation 1, 
while it reaches only x ~ —50 and v x ~ — 2t> i2 in simulation 2. This can be attributed 
to the better representation of the high-energy tail of the Maxwellian in simulation 1. 

The bulk electron distributions in both simulations agree well for x < 100. The 
interaction of the confined electrons of plasma 1 with the expanding protons is thus 
reproduced well by both simulations. We find a beam of electrons with x > 100 and 
v x w — 3t> t i in Fig. [5|c), which is accelerated by the double layer to — 4t> tl in the interval 
—100 < x < 100. This beam originates from the second boundary between the dense 
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Figure 6. (Colour online) The 10-logarithmic phase space distributions of the 
snapshots Si (a,b) and S 2 (c,d) in units of CPs: A shock is developing in the proton 
distribution (a) at x « 300. The electrons are distributed symmetrically around v x = 
in (b) and their density value jumps at x » 300. The proton shock in (c) and the 
electron density jump in (d) have propagated to x ~ 600 at the mean speed « v t 2- 



and the tenuous plasma at a; = L$/2 in simulation 1. It is thus an artifact of our periodic 
boundary conditions. Its density is three orders of magnitude below the maximum one 
and it does thus not carry significant energy. This tenuous beam does not show any 
phase space structuring, which would be a consequence of instabilities, and it has thus 
not interacted with the bulk plasma. Its only consequence is to provide a weak current 
that should not modify the double layer. This fast beam is absent in Fig. 0(d), because 
the electrons could not cross the distance Ll/2 in simulation 2 during the time tg- 

The electron distributions for x > 100 and v x > computed by both simulations 
differ substantially. The electrons form phase space vortices in simulation 1, while the 
electrons in simulation 2 form a diffuse beam with some phase space structures, e.g. 
at x ~ 300. Phase space vortices are a consequence of an electrostatic two-stream 
instability, which must have developed between the leaked electrons of plasma 1 and 
the electrons of plasma 2. Only the electrons of plasma 1 with v > 3v t ± can overcome 
the double layer potential. These leaked electrons form a smooth beam in simulation 1 
that can interact resonantly with the electrons of plasma 2 to form well-defined phase 
space vortices. The statistical representation of the leaking electrons in simulation 2 
provides a minimum density that exceeds the density of these vortices. 

4. Simulation 2: Long term evolution 

We examine the plasma at three times. The snapshot Si corresponds to the time 
t = 8000, 5*2 to t = 16000 and S3 to t = 25500. The plasma phase space distributions 
for Si and S2 are displayed by Fig. El The proton distribution is still qualitatively 
similar to that at t — 300 in Fig. HI The phase space boundary between the protons 
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of plasma 1 and 2 has been tilted further by the proton streaming. The key difference 
between the Figs. H] and [6] is found, where the proton distribution of the plasma 1 
merges with that of the plasma 2. This collision boundary is located at x ~ 300 for 
Si and at i « 600 for S 2 , which evidences an approximately constant speed of this 
intersection point. The propagation speed is ~ Vt 2 . The protons directly behind this 
collision boundary, e.g. in 450 < x < 550 for S 2 , do not show a velocity shear. Their 
mean speed and velocity spread is spatially uniform in this interval, evidencing the 
downstream region of a shock. The upstream proton distribution with x > 600 for S 2 
resembles, however, only qualitatively that of an electrostatic shock [11] . That consists 
of the incoming plasma and the shock-reflected ion beam. The density of the beam with 
v x ~ Av t2 exceeds that of the plasma 2 in the same interval and its mean speed exceeds 
the v s ~ Vt2 of the shock by the factor 4. A shock-reflected ion beam would move at 
twice the shock speed and its density would typically be less than that of the upstream 
plasma, which the shock reflects. The linear increase of the proton beam velocity with 
increasing x is reminiscent of the plasma expansion into a vacuum [20], but it is here a 
consequence of the shear introduced by the proton thermal spread. 

The electron distribution at t = t$ in Fig. E(d) could be subdivided into the cool 
electrons of plasma 2 and the leaked hot electrons of plasma 1, while the electrons in 
the interval x > 750 have a symmetric velocity distribution in Fig. MJd) that does not 
permit such a distinction. The electron temperature gradient has also been eroded. 
The electron phase space density decreases by an order of magnitude as we go from 
v x = to v x ~ 2vn at x ~ and at x » 2000 in Fig. 0(d) and the thermal spread 
is thus comparable at both locations. We attribute this temperature equilibration to 
electrostatic instabilities, which were driven by the electron beam that leaked through 
the boundary at x — 0, and also to the electron scattering by the simulation noise. 
The noise amplitude is significant in the interval x > due to the comparatively low 
statistical representation of the plasma, in particular that of the hot leaked electrons. 

The electron density jumps at both times in Fig. [6] at the positions, where the 
protons of plasma 1 and 2 intersect. The electron distribution for S 2 furthermore shows 
a spatially uniform distribution in 450 < x < 550, as the protons do. The electrons have 
thermalized and any remaining free energy would be negligible compared to that of the 
protons. The electron density merely follows that of the protons to conserve the plasma 
quasi-neutrality. This electron distribution thus differs from the similarly looking one, 
which has been computed at late times in Ref. [21] . There the electrons changed their 
velocity distribution in response to the energy lost to the protons. 

The time 10t$ corresponding to Si and the box length Ll = IOL5 imply, that we 
should see some electrons emanated by the plasma boundary at x = Ljj/2 as in Fig. 
[5j Only the electrons with v < —2.1v t i would be fast enough to cross the interval 
< x < Ll/2 occupied by plasma 2 during the time lOts- These electrons correspond 
to the few fast electrons in Fig. Mjo) with x > and v < 0. An increased number 
of fast electrons moving in the negative x-direction is visible at the snapshot S 2 - The 
electrons emanated from the plasma boundary at x = Ll/2 reach now the boundary 
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(b) Position (d) Position 

Figure 7. (Colour online) The 10-logarithmic phase space density for S3 in units of 
CPs: (a) displays the proton distribution and (b) the electron distribution. The shock 
is located at x ks 900 and phase space holes develop in the proton (c) and electron (d) 
distribution at 1600 < x < 1700. A new shock is growing at x ks 1700 in (c). 



at x = in significant numbers. The diffuse phase space distribution of these electrons 
implies, however, that they do not carry with them enough free energy that could result 
in instabilities that drive strong electrostatic fields. 

The shock structure and the density jump in the electron distribution has 
propagated to x ~ 900 for S3 and the proton beam ahead of the shock has started 
to thermalize by its interaction with the upstream plasma, as it is evidenced by the Fig. 
El An electron phase space hole doublett and proton phase space structures are visible. 
These structures have grown out of the phase space oscillation of the proton beams and 
the electron phase space hole at x ~ 1250 in Fig. [6](c). The proton distribution in 
Fig. [3(c) in x > 1700 reveals that a second shock is forming, which will thermalize the 
dense and fast beam of protons that expands out of the plasma 1 into the plasma 2. 
The spatially uniform electron distribution outside the interval occupied by the electron 
phase space holes changes only its thermal spread and density along x and could be 
approximated by a Boltzmann-distribution. The electrons are not accelerated to high 
energies neither by the shocks nor by the other phase space structures. 

The expansion of the protons of plasma 1 in simulation 2 is captured by the movie 
3. The colour scale corresponds to the 10-logarithmic number of CPs. The movie 3 
evidences the formation of the shock and of its downstream region and it displays how 
the proton phase space hole and, subsequently, the secondary shock develop. The mean 
velocity of the upstream protons is modulated along x, which is probably a result of the 
same wave fields that thermalized the electrons. 

The proton distribution at x ~ changes in time primarily due to the free motion 
of a proton % with the speed v xil which is displaced clS X{ — V x it. The phase space 
boundary between the plasma 1 and 2 is thus increasingly sheared. Further acceleration 
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Scaled position 

Figure 8. The proton densities, normalized to rih, as a function of the scaled position 
xti/tj, where tj corresponds to the snapshot Sj. The curves match, except within 
the downstream region of the shock at 200 < xti/tj < 400 that is characterized by a 
constant density. The density doubles by the shock compression at xti/tj w 350. 




(b) Position (d) Position 

Figure 9. (Colour online) The 10-logarithmic electron phase space distributions are 
shown in (a) for Si and (c) for S3. The electron density (dashed curves) and the 
electrostatic field (solid curves) are displayed in (b) for Si and in (d) for S3. The 
densities are integrated and the electric fields averaged over 5 grid cells. 

mechanisms are the drag of the protons by the thermally expanding electrons and 
the shock formation. Figure M assesses their relative importance. The plasma density 
distribution should be invariant if the protons expand freely and if we scale the position 
oc x/t. This is indeed the case and the proton density distributions for Si,S% and S3 
match if we use the scaled positions, except at the shock and within its downstream 
region. The electron densities (not shown) closely follow those of the protons. 

Figure [9] compares the electrostatic field with the electron distributions for the 
snapshots Si and S3. An electric field peak at x ~ 400 coincides with the shock in the 
snapshot S\. The peak E x « 0.04 and it confines the electrons to the left of the density 
jump by accelerating them into the negative x direction. The electric field can be scaled 
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to physical units with tim — 10 15 cm -3 and v t \ = 1.325 x 10 7 m/s to give ~ 5 x 10 6 V/m. 
The electric field, which has been measured close to the shock in Ref. [8], is < 2 x 10 7 
V/m. The plasma density in the region, where the shock develops in the experiment, 
may be higher than 10 15 cm" 3 . The electric field amplitudes associated with the shock 
are thus comparable. The noise levels in PIC simulations are typically higher than in 
a physical plasma, explaining the strength of the evenly spread noise in the simulation 
box, which is not observed to the same extent in the experiment. The electric field at 
the shock at x ~ 10 3 is at noise levels for S3, while the phase space holes at x ~ 1700 
give an electric field, which is exceeding that sustained by the shock for Si. 

5. Discussion 

We have investigated the thermal expansion of a hot dense plasma into a cooler tenuous 
plasma. The thermal pressure of the hot plasma exceeded that of the cool plasma by 
the factor 10 4 . Our study has been motivated by the laser-plasma experiment in Ref. 
[E], which examined the expansion of a hot and dense plasma into a tenuous ambient 
medium. Our initial conditions and the ID geometry are, however, idealized and the 
simulation results can thus not be compared quantitatively to the experimental ones. 
The aim of our work has been to better understand the qualitative effects of the ambient 
medium on the plasma expansion. We have for this purpose compared our results with 
some of those in the related study in Ref. [21], that considered the plasma expansion 
into a vacuum. There, the electron temperature exceeded that of the protons by a factor 
10 3 , while we consider here the same temperature of electrons and protons. 

Our results are summarized as follows. An electric field grows almost instantly at 
the boundary between both plasmas, because the ion expansion of the hot plasma is 
slower than the electron expansion. The electric field forms irrespective of the ambient 
medium. It accelerates only the ions, if the plasma expands into a vacuum and it has 
a cusp in its spatial profile. The acceleration of the electrons of the ambient medium 
triggers in our simulation the formation of a double layer [22] with a smooth electric 
field profile. This double layer redistributes the momentum between the individual 
plasma species [25J . A tenuous hot beam of electrons streams from the hot plasma into 
the cool plasma, while all the electrons of the cool plasma are dragged into the hot 
plasma. These beams thermalize through electrostatic two-stream instabilities, which 
equilibrate the electron temperatures of both plasmas on electron time scales. This rapid 
thermalization will cancel any significant proton acceleration by hot electrons already 
at the relatively low density of the ambient medium we have used. Proton acceleration 
is, however, still possible because a thermal pressure gradient is provided by the density 
jump. Most electrons merely follow after their thermalization the motion of the protons 
to conserve the quasi-neutrality of the plasma. They maintain their Maxwellian velocity 
distribution, which would not be the case for an expansion into a vacuum [2T] . 

The protons at the front of the hot plasma are accelerated by the electric field of 
the double layer to about 5.5 times the proton thermal speed, while the Maxwellian 
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distribution is represented up to 3-4 times the proton thermal speed. The expansion of 
the protons from the hot into the cool plasma is dominated by the free streaming of the 
fastest protons (diffusion). The effects of the ambient medium on the proton expansion 
are initially negligible. Eventually the interaction of the expanding and the ambient 
plasma results in the formation of shocks. We have observed one shock at the position, 
where the protons of both plasmas merge. This shock did not result in the acceleration 
of electrons or in the modification of their phase space distribution. 

The protons of the hot plasma expand farther than the position of this shock and 
they can interact with the protons of the cool plasma through ion beam instabilities. 
The interval, in which the protons of both plasmas co-exist, resembles qualitatively the 
upstream region of an electrostatic shock [IT] . However, the density and the speed of the 
beam of expanding protons of the hot plasma are both higher than what we expect for 
a shock-reflected ion beam. We have observed in the simulation the growth of a phase 
space structure in the upstream proton distribution that gave rise to an electron phase 
space hole. The proton structure evolved into a second shock ahead of the primary one. 
The presence of multiple shocks has been observed experimentally [26], although there 
the second shock was radiation-driven and not beam-driven. 
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